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Abstract 

We present an algorithm which calculates groundstates of Ising spin 
glasses approximately. It works by randomly selecting clusters of spins 
which exhibit no frustrations. The spins which were not selected, con- 
tribute to the local fields of the selected spins. For the spin-cluster a 
groundstate is exactly calaculated by using graphtheoretical methods. 
The other spins remain unchanged. This procedure is repeated many 
times resulting in a state with low energy. The total time complexity of 
this scheme is approximately cubic. We estimate that the groundstate 
energy density of the infinite system for the ± J model is —1.400 ± 0.005 
(2d) and —1.766 ± 0.002 (3d). The distribution of overlaps for selected 
systems is calculated in order to characterize the algorithm. 



The combination of frustration and randomness makes it difficult to find ground- 
states of spin glasses M using numerical simulations. In the past years many 
methods @ have been proposed including the multicanonical ensemble ||, ge- 
netic algorithms a scheme, which uses storing of spin configurations [||], and 
an exact algorithm exhibiting exponential timecomplexity 0. In this paper 
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we present a method which combines clustering of spins with exact calculation 
of groundstates of polynomial solvable problems in order to calculate approxi- 
mately groundstates of Edwards- Anderson (EA) spin glasses. 
The Hamiltonion of the EA model is given by 



H = - JijOiVj - BiOi (1) 

<ij> i 

where the sum goes over nearest neighbours of spins Oi = ±1. The exchange 
interactions Jy are selected at random according to a probability distribution. 
The external field Bi can be site-dependend. 

For our computer experiments we used simple cubic lattices of N = L d spins 
(L = linear lenght, d = dimension) with periodic boundary conditions in no 
external field (Bi = 0). The interactions are Jy = ±1 with equal probability 
using the constraint J2<ij> Jij = 0. 

The algorithm for the approximation of groundstates in frustrated Ising systems 
works by taking a spin configuration and calculating another one, which has a 
lower or equal energy. This procedure is iterated many times. 
The idea of the scheme for lowering the energy is to choose a cluster of spins, 
which exhibits no frustrations. The interaction of the cluster with the spins at its 
boundary is included into the local fields of the cluster-spins. The groundstate 
of the cluster is exactly calculated by using concepts of graph theory [|7|, g, [|: 
An equivalent network is constructed |l(J , the maximum flow is calculated with 
the Ford-Fulkerson algorithm [IT| and a minimum cut is constructed [O. The 
new configuration consists of the unchanged spins, which are not in the cluster, 
and the groundstate of the cluster. By definition this procedure cannot increase 
the energy. 

For describing the spin-cluster the Hamiltonian can be rewritten in the form 



He = - £ JvUtrfof - £ + C (2) 

<ij> i 

The values of t\ = 0, ±1 describe the cluster and are used to handle antiferro- 
magnetic interactions. If spin (jj =: tjCr| does not belong to the cluster: ^ = 0. 
If ti and tj 7^ their signs have to be choosen so that J?- := Jijitfj > 0, be- 
cause only Hamiltonians with positive exchange interactions can be transformed 
into an equivalent network. The local fields B c { include the external field Bi and 
the interactions with the neighbours of spin (jj, which are not in the cluster: 

Bt = Bi + Y. J iA±-M)°i (3) 

<j> 

The constant C summarizes the interactions between the non-cluster spins and 
can be dropped. The construction of the unfrustrated cluster works in the 
following way: a spin is randomly selected as seed. Iteratively neighbouring 
spins are added, if no frustration occours. If a cluster cannot be extended, one 
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more seed-spin is selected. The following algorithm contains the details. The 
local variables Si are used to indicate the spins, which are already detected as 
neighbours of the cluster. The set A contains the neighbours of spin a i: which 
are already in a cluster, whereas B contains the neighbours, which are added 
to the boundary of the cluster. 

algorithm create_cluster({ J y -}) 
begin 

for all i do 
begin 

initialise U <— 0; 
initialise Si <— 0; 

end; 

while there are unmarked spins with Si = do 
begin 

select one index i from the spins with Si = 0; 
push i on an empty Stack S; 
k <- l; 

while S is not empty do 
begin 

pop one randomly choosen index % from S; 
A <— {j\j is neighbour of % and tj ^ 0}; 
if A = then t t <- +1; 

else if Vj G A : Jijtj has the same sign a then 

begin 

ti < <y, 

B <— {j\j is neighbour of i and Sj = 0}; 
for all j e B do 

begin 

push j on 5; 

Sj - l; 

end; 

end; 

else ti <— 0; 

end; 

end; 

return ({*»}); 

end; 

A run for calculating a groundstate consists of choosing randomly an initial 
configuration or choosing all spins pointing up, and calculating new configura- 
tions until the energy can not further be lowered. For this in our experiments 
we found as a good criterium that the energy did not change for the last n g 
steps. We used n g = N/2 (2d, 3d) as rule of thumb, because in some very long 
runs we observed, that the longest period between two jumps in energy never 
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exceeded this value and scaled almost linear with system size. 
Because the minimum energy reached in a run depends on the starting config- 
uration and on the clusters, which were constructed, we performed several runs 
for each realization of the random variables Jy. We used that one, which ended 
on the lowest energy level. We found 3 runs per realization sufficent, because by 
increasing the number of runs, the average groundstate energy was only lowered 
about 0.01 percent. 

Because in each step the spin-cluster is randomly built and by using the algo- 



rithm |L2[ all of the degenerated cluster-groundstates have a positive probability 
to be calculated, usually in each step a new spin configuration is constructed, 
even if the energy remains constant. So it is possible to explore large areas of 
the configuration space. 

In figure [l] the average energy density ez,(£) = Ei{t)/N of 96 realizations is 
shown as a function of the step number for lattice sizes L = 4,6, 16 of the 2d 
system. One can see, that in the beginning the algorithm approaches very fast 
low values. Later the decrease in energy is very small. 

It is possible to use other procedures in selecting the spins of the cluster. We 
tried three other methods (we call the first presented algorithm method A): 

method B Not only neighbours of cluster-spins were added to the stack S, 
but neighbours of all spins, which were tested, if they can be 
added to the cluster. 

method C The spins which were tested, if they can be added to the cluster 
were totaly selected at random from the spins, which were yet 
untested. 

method D Methods A and B are applied alternately. 

We tested all four methods by calculating groundstates for eight realizations of 
16 2 systems (3 runs each). In table |l| the average groundstate energy density 
e° is displayed. Also we constructed for the 2d and 3d case 100 respectively 10 
spin-clusters for 100 realizations of system sizes L — 4,6, ... , 20. The average 
size fractions n c /N of the spin-clusters are shown in the last two rows. 



method 


A 


B 


C 


D 


e u 


1.425 


1.425 


1.415 


1.427 


n c /N (2d) 


0.708 


0.696 


0.644 




n c /N (3d) 


0.578 


0.568 


0.541 





Table 1: Comparison of four methods used to construct spin-clusters 



We observed, that for some realizations method A and for others method B 
gave lower energies, so it is plausible, that the combination of both methods 
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gave the best results. So we used method D for all further experiments. Maybe 
better methods for constructing the cluster would result in a faster convergence 
and/or lower energies. 

To estimate the groundstate energies of the 2d and 3d ±J spin glass we per- 
formed calculations for lattice sizes L — 4, ... ,30 (2d) and L — 4, ... ,14 (3d) 
with 96 (4 2 , . . . , 20 2 , 4 3 , . . . , 10 3 ), respectively 32 (30 2 , 12 3 , 14 3 ) realizations of 
the random variables Jy. 

In figure |2| the groundstate energy density e°(L) of the 2d systems is plotted 
as function of system size. To estimate the groundstate energy density of the 
infinite system we performed a finite-size scaling analysis. A fit to the function 
fl(N) = + c/N m results in e° = -1.400 ± 0.005. This value is consistent 
with earlier results from a genetic algorithm —1.400 ± 0.005 M, multicanonical 



simulations -1.394±0.007 @ pure MC -1.407±0.008 [13 and transfer matrix 



calculations -1.4024 ± 0.0012 HI 



The results for the 3d systems are shown in figure [| The fit gives a groundstate 
energy density of e° = — 1.766±0.002. The results of other authors are 1.7863± 
0.0028 and 1.765 ± 0.01 f| 

As seen in figure |l] the time for approaching the groundstate increases rapidly 
with the system size. The average number of steps needed to calculate a ground- 
state are shown in figure § as a function of the number of spins N. For the 2d 
system (lower curve) a fit to a function t(N) = a + bN c results in a timecom- 
plexity of O(n 117±006 ). The result for the 3d system is similiar: O(n L27±a08 ). 
Because the time to calculate one groundstate of the cluster increases quadrat- 



icly with the number of spins |IIJ and the size of the cluster is linear in the 
system size, this totally results in an approximate cubic timecomplexity of the 
algorithm. Performing one run for example for a 10 3 system consisting of 1000 
configurations needed about 10 hours on a 80 Mhz PowerPC processor system. 
In order to characterize the way the algorithm approaches the groundstates we 
selected one realization of a 2d (L = 16) system and used configurations of 
different energy densities e s G [—1.414, 1.0] as starting-configurations. For each 
configuration 8 different runs for constructing groundstates were performed. For 
each run, which resulted in the groundstate energy density value of e° = —1.414 
the last 100 configurations were stored. The overlap 

N 

?:=£°f*f ( 4 ) 



was calculated between all pairs (a,/3) of groundstates belonging to a starting- 
configuration . We obtained the corresponding probability distributions Pl{q) 
by counting the numbers of overlaps within intervals of length Aq = 0.05. The 
result is displayed in figure ^]for the starting energies e s = —1.414, —1.07, 0.117. 
For these three staring energies 800, 300 respectively 400 of the resulting states 
exhibited the lowest energy density. One can see, that the algorithm approaches 
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local minima of the energy landscape. The higher the starting energy the more 
minima are reachable. If a system is caught in a local minimum, only a restricted 
area of configuration space can be explored by the algorithm, although in each 
step a different configuration is generated which is indicated by -Pl(I) = 0. 
But for other realizations we got also different distributions. There are systems 
where Pl{q) is always narrow or always broad, independent of the starting 
energy. The accessible area is usually smaller for lower groundstate energies 
than for higher lying groundstates. 

In conclusion, in this letter we have presented an algorithm for the approxima- 
tion of groundstates for Ising spin glasses. Its time complexity is approximately 
cubic in the number of spins. Because by each step a large number of spins is 
allowed to flip, the algorithm approaches low energies very rapidly. 
So it should be interesting to try combinations with other algorithms like ge- 
netic algorithms or Monte Carlo simulations, or to use it for other Ising type 
optimization problems. 

Although we applied it to the ± J model it is also possible to treat more general 
forms of lattices, interactions and their probability distributions. 
It should possible to get results for B ^ in a resonable amount of time in 
order to characterize better the phase diagramm along the T = axis. 
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Figure Captions 



Figure [I] Average energy density e^(t) = E L (t)/N of 2d ±J spin glass at different 
iterations t. Shown are 3 sizes L = 4, 6, f 6. 

Figure |2] Average groundstate energy density e°(L) = E°(L)/N of 2d ± J spin glass 
function of size L. 

Figure |3] Average groundstate energy density e°(L) = E°(L)/N of 3d ±J spin glass 
as a function of size L. 

Figure ^ Average number of steps needed to reach a groundstate as function of 
system size N for 2d system (L = 4, . . . , 20). 

Figure |5] Distribution of overlaps for one single realization (2d, L = 16, e° = 

— 1.414). For different starting configurations with energies e s = 0.117, 

— 1.07, —1.414 with 8 independent runs 100 groundstates were calculated. 
The states with energy density e° = —1.414 were included in the calcula- 
tion of the distribution, i.e. 400, 300, 800 states for the three values of e s . 
The lines are guides for the eyes only. 
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